Model Diagnostic Comparisons

Overview

Diagnostics from initial models demonstrated moderate issues with kurtosis, nonlinearity, and heteroskedasticity. The goal of this comparison is to explore whether other modeling options improve these issues. The first focal species is vermillion rockfish (Sebastes miniatus), a species found statewide that is generally associated with hard bottom habitat. This was chosen because it has two int

The original approach was:

Generalized linear mixed effects models with all numeric variables scaled (including dep. variable, first log+1 transformed then scaled to mean = 0 and sd = 1). The variables were scaled to accommodate differences in scales between MPA age and the habitat variables, and the magnitude between the various buffer scales now that multiple scales could be included in the same models (e.g. hard bottom at 500m and kelp at 25m, varying from 10s of sq meters to millions of sq meters).

The top model for this species was:

Log+1-Transformed Biomass ~ Hard Bottom 500m * Site Type + Kelp Annual 100m * Site Type + Site Type * MPA Age

-   Does log-transforming biomass make a difference in model fit? Is there a difference between using no transformation, log+1, vs. log+c, where c is a very small constant (10x smaller than the smallest nonzero vale)

-   Does scaling the variables (mean = 0 and sd = 1) make a difference in conclusions or model fit?

-   Does adding a 3-way interaction with MPA age make a difference in model fit?

Build

# Define constant 10x smaller than minimum nonzero biomass
const <- min(data_kelp$kg_per_m2[data_kelp$kg_per_m2 > 0])/10

data_sp <- data_kelp  %>% 
  mutate(
    # Add small constant 
    kg_per_m2_c = kg_per_m2 + const,
    # Log-transform to variable with small constant
    log_kg_per_m2_c = log(kg_per_m2_c))

# Select term columns to join with residuals df
data_terms <- data_sp %>% 
  dplyr::select(obs_id, age_at_survey, site_type, 
                all_of(grep("^(hard|soft|kelp|depth)", names(.), value = TRUE)))

data_sp_scaled <- data_sp %>% 
  mutate(across(where(is.numeric), ~ scale(.) %>% as.vector()))

Untransformed

B ~ H*ST + K*ST + ST*A (Not Scaled)

  • Removed year as it had zero variance (singular fit)
  • Significant: H*ST, K*ST, ST*A, ST
glm1 <- lmer(kg_per_m2 ~ hard_bottom_500 * site_type + kelp_annual_100 * site_type + site_type * age_at_survey +
               (1|bioregion) + (1|affiliated_mpa), data = data_sp, REML = FALSE)
Warning: Some predictor variables are on very different scales: consider
rescaling
Warning: Some predictor variables are on very different scales: consider
rescaling
attr(glm1, "name") <- "GLMM: B ~ H*ST + K*ST + ST*A (Not Scaled)"

summary(glm1)
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: 
kg_per_m2 ~ hard_bottom_500 * site_type + kelp_annual_100 * site_type +  
    site_type * age_at_survey + (1 | bioregion) + (1 | affiliated_mpa)
   Data: data_sp

     AIC      BIC   logLik deviance df.resid 
-15536.1 -15475.9   7779.0 -15558.1     1737 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6605 -0.3363 -0.0311  0.1542 14.3967 

Random effects:
 Groups         Name        Variance  Std.Dev. 
 affiliated_mpa (Intercept) 9.310e-07 0.0009649
 bioregion      (Intercept) 1.968e-06 0.0014028
 Residual                   7.744e-06 0.0027829
Number of obs: 1748, groups:  affiliated_mpa, 23; bioregion, 3

Fixed effects:
                               Estimate Std. Error         df t value Pr(>|t|)
(Intercept)                   1.975e-03  9.058e-04  3.631e+00   2.180 0.101628
hard_bottom_500              -9.557e-10  5.811e-10  1.630e+03  -1.645 0.100214
site_typeMPA                 -3.139e-03  4.155e-04  1.721e+03  -7.555 6.75e-14
kelp_annual_100              -3.892e-08  2.924e-08  1.745e+03  -1.331 0.183258
age_at_survey                 1.928e-05  1.785e-05  1.743e+03   1.080 0.280182
hard_bottom_500:site_typeMPA  1.009e-08  1.080e-09  1.572e+03   9.342  < 2e-16
site_typeMPA:kelp_annual_100  1.529e-07  3.639e-08  1.745e+03   4.203 2.77e-05
site_typeMPA:age_at_survey    9.706e-05  2.567e-05  1.734e+03   3.781 0.000161
                                
(Intercept)                     
hard_bottom_500                 
site_typeMPA                 ***
kelp_annual_100                 
age_at_survey                   
hard_bottom_500:site_typeMPA ***
site_typeMPA:kelp_annual_100 ***
site_typeMPA:age_at_survey   ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) hr__500 st_MPA k__100 ag_t_s h__500: s_MPA:__1
hrd_btt_500 -0.234                                               
site_typMPA -0.146  0.297                                        
klp_nnl_100 -0.124  0.130   0.197                                
age_at_srvy -0.185  0.061   0.361  0.246                         
h__500:_MPA  0.041 -0.309  -0.740 -0.003  0.033                  
s_MPA:__100  0.099 -0.124  -0.255 -0.732 -0.200 -0.037           
st_tyMPA:__  0.139 -0.074  -0.521 -0.161 -0.681 -0.050   0.215   
fit warnings:
Some predictor variables are on very different scales: consider rescaling

B ~ H*ST + K*ST + ST*A (Scaled)

  • Removed year as it had zero variance (singular fit)
  • Significant: H*ST, K*ST, ST*A, ST
  • Same as nonscaled version (#1)
glm2 <- lmer(kg_per_m2 ~ hard_bottom_500 * site_type + kelp_annual_100 * site_type + site_type * age_at_survey +
               (1|bioregion) + (1|affiliated_mpa), data = data_sp_scaled, REML = FALSE)

attr(glm2, "name") <- "GLMM: B ~ H*ST + K*ST + ST*A (Scaled)"

summary(glm2)
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: 
kg_per_m2 ~ hard_bottom_500 * site_type + kelp_annual_100 * site_type +  
    site_type * age_at_survey + (1 | bioregion) + (1 | affiliated_mpa)
   Data: data_sp_scaled

     AIC      BIC   logLik deviance df.resid 
  4306.8   4366.9  -2142.4   4284.8     1737 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6605 -0.3363 -0.0311  0.1542 14.3967 

Random effects:
 Groups         Name        Variance Std.Dev.
 affiliated_mpa (Intercept) 0.07924  0.2815  
 bioregion      (Intercept) 0.16748  0.4092  
 Residual                   0.65916  0.8119  
Number of obs: 1748, groups:  affiliated_mpa, 23; bioregion, 3

Fixed effects:
                               Estimate Std. Error         df t value Pr(>|t|)
(Intercept)                     0.19558    0.25263    3.04429   0.774 0.494481
hard_bottom_500                -0.04452    0.02707 1630.32168  -1.645 0.100214
site_typeMPA                    0.28866    0.04080 1742.29575   7.075 2.16e-12
kelp_annual_100                -0.04514    0.03390 1744.51016  -1.331 0.183258
age_at_survey                   0.03029    0.02804 1743.12071   1.080 0.280182
hard_bottom_500:site_typeMPA    0.47000    0.05031 1571.97148   9.342  < 2e-16
site_typeMPA:kelp_annual_100    0.17734    0.04220 1745.12473   4.203 2.77e-05
site_typeMPA:age_at_survey      0.15247    0.04032 1734.11137   3.781 0.000161
                                
(Intercept)                     
hard_bottom_500                 
site_typeMPA                 ***
kelp_annual_100                 
age_at_survey                   
hard_bottom_500:site_typeMPA ***
site_typeMPA:kelp_annual_100 ***
site_typeMPA:age_at_survey   ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) hr__500 st_MPA k__100 ag_t_s h__500: s_MPA:__1
hrd_btt_500 -0.024                                               
site_typMPA -0.080 -0.007                                        
klp_nnl_100  0.003  0.130  -0.030                                
age_at_srvy  0.007  0.061   0.008  0.246                         
h__500:_MPA -0.013 -0.309   0.026 -0.003  0.033                  
s_MPA:__100 -0.003 -0.124  -0.021 -0.732 -0.200 -0.037           
st_tyMPA:__  0.003 -0.074  -0.024 -0.161 -0.681 -0.050   0.215   

B ~ H*ST*A + K*ST*A (Not Scaled)

  • Removed year as it had zero variance (singular fit)
  • Significant: H*ST*A, K*ST, ST*A
  • Almost significant: H, ST
glm3 <- lmer(kg_per_m2 ~ hard_bottom_500 * site_type * age_at_survey + kelp_annual_100 * site_type * age_at_survey + 
               (1|bioregion) + (1|affiliated_mpa), data = data_sp, REML = FALSE)
Warning: Some predictor variables are on very different scales: consider
rescaling
Warning: Some predictor variables are on very different scales: consider
rescaling
attr(glm3, "name") <- "GLMM: B ~ H*ST*A + K*ST*A (Not Scaled)"

summary(glm3)
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: 
kg_per_m2 ~ hard_bottom_500 * site_type * age_at_survey + kelp_annual_100 *  
    site_type * age_at_survey + (1 | bioregion) + (1 | affiliated_mpa)
   Data: data_sp

     AIC      BIC   logLik deviance df.resid 
-15611.3 -15529.3   7820.6 -15641.3     1733 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.9976 -0.3236 -0.0405  0.1645 14.7939 

Random effects:
 Groups         Name        Variance  Std.Dev. 
 affiliated_mpa (Intercept) 9.166e-07 0.0009574
 bioregion      (Intercept) 2.114e-06 0.0014539
 Residual                   7.381e-06 0.0027167
Number of obs: 1748, groups:  affiliated_mpa, 23; bioregion, 3

Fixed effects:
                                             Estimate Std. Error         df
(Intercept)                                 2.155e-03  9.572e-04  3.952e+00
hard_bottom_500                            -1.674e-09  9.524e-10  1.744e+03
site_typeMPA                                4.896e-04  6.226e-04  1.739e+03
age_at_survey                              -2.231e-06  3.181e-05  1.740e+03
kelp_annual_100                            -2.060e-08  4.054e-08  1.746e+03
hard_bottom_500:site_typeMPA               -1.416e-09  1.875e-09  1.717e+03
hard_bottom_500:age_at_survey               8.609e-11  9.125e-11  1.736e+03
site_typeMPA:age_at_survey                 -3.415e-04  6.185e-05  1.737e+03
site_typeMPA:kelp_annual_100                1.342e-07  5.380e-08  1.737e+03
age_at_survey:kelp_annual_100              -3.951e-09  5.866e-09  1.744e+03
hard_bottom_500:site_typeMPA:age_at_survey  1.380e-09  1.870e-10  1.739e+03
site_typeMPA:age_at_survey:kelp_annual_100  2.683e-09  7.479e-09  1.731e+03
                                           t value Pr(>|t|)    
(Intercept)                                  2.251   0.0884 .  
hard_bottom_500                             -1.758   0.0790 .  
site_typeMPA                                 0.786   0.4318    
age_at_survey                               -0.070   0.9441    
kelp_annual_100                             -0.508   0.6114    
hard_bottom_500:site_typeMPA                -0.755   0.4504    
hard_bottom_500:age_at_survey                0.943   0.3456    
site_typeMPA:age_at_survey                  -5.522 3.86e-08 ***
site_typeMPA:kelp_annual_100                 2.494   0.0127 *  
age_at_survey:kelp_annual_100               -0.674   0.5007    
hard_bottom_500:site_typeMPA:age_at_survey   7.379 2.45e-13 ***
site_typeMPA:age_at_survey:kelp_annual_100   0.359   0.7199    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
             (Intr) hr__500 st_MPA ag_t_s k__100 hr__500:_MPA h__500:__
hrd_btt_500  -0.316                                                    
site_typMPA  -0.175  0.398                                             
age_at_srvy  -0.292  0.669   0.429                                     
klp_nnl_100  -0.130  0.118   0.185  0.309                              
hr__500:_MPA  0.119 -0.431  -0.889 -0.318 -0.053                       
hrd__500:__   0.234 -0.803  -0.353 -0.813 -0.094  0.405                
st_tyMPA:__   0.154 -0.350  -0.830 -0.501 -0.155  0.736        0.413   
s_MPA:__100   0.094 -0.088  -0.179 -0.229 -0.714 -0.040        0.068   
ag_t_:__100   0.054 -0.035  -0.108 -0.248 -0.708  0.043        0.065   
h__500:_MPA: -0.116  0.391   0.740  0.387  0.050 -0.826       -0.481   
s_MPA:__:__  -0.037  0.018   0.069  0.187  0.536  0.058       -0.050   
             st_MPA:__ s_MPA:__1 a__:__ h__500:_MPA:
hrd_btt_500                                         
site_typMPA                                         
age_at_srvy                                         
klp_nnl_100                                         
hr__500:_MPA                                        
hrd__500:__                                         
st_tyMPA:__                                         
s_MPA:__100   0.122                                 
ag_t_:__100   0.123     0.530                       
h__500:_MPA: -0.897     0.065    -0.033             
s_MPA:__:__  -0.060    -0.748    -0.761 -0.113      
fit warnings:
Some predictor variables are on very different scales: consider rescaling

B ~ H*ST*A + K*ST*A (Scaled)

  • Removed year as it had zero variance (singular fit)
  • Significant: H*ST*A, H*ST, K*ST, ST*A, ST
  • Almost significant: H
  • Similar to non-scaled version except site type is now highly significant on its own
glm4 <- lmer(kg_per_m2 ~ hard_bottom_500 * site_type * age_at_survey + kelp_annual_100 * site_type * age_at_survey + 
               (1|bioregion) + (1|affiliated_mpa), data = data_sp_scaled, REML = FALSE)

attr(glm4, "name") <- "GLMM: B ~ H*ST*A + K*ST*A (Scaled)"

summary(glm4)
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: 
kg_per_m2 ~ hard_bottom_500 * site_type * age_at_survey + kelp_annual_100 *  
    site_type * age_at_survey + (1 | bioregion) + (1 | affiliated_mpa)
   Data: data_sp_scaled

     AIC      BIC   logLik deviance df.resid 
  4231.6   4313.6  -2100.8   4201.6     1733 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.9976 -0.3236 -0.0405  0.1645 14.7939 

Random effects:
 Groups         Name        Variance Std.Dev.
 affiliated_mpa (Intercept) 0.07802  0.2793  
 bioregion      (Intercept) 0.17991  0.4242  
 Residual                   0.62820  0.7926  
Number of obs: 1748, groups:  affiliated_mpa, 23; bioregion, 3

Fixed effects:
                                             Estimate Std. Error         df
(Intercept)                                   0.18751    0.26058    3.00887
hard_bottom_500                              -0.04453    0.02646 1635.65954
site_typeMPA                                  0.28244    0.04112 1741.50636
age_at_survey                                 0.02530    0.02866 1743.68488
kelp_annual_100                              -0.06210    0.04063 1743.73621
hard_bottom_500:site_typeMPA                  0.47020    0.04928 1585.04419
hard_bottom_500:age_at_survey                 0.02159    0.02289 1736.27602
site_typeMPA:age_at_survey                    0.12217    0.04042 1734.05950
site_typeMPA:kelp_annual_100                  0.18151    0.04876 1743.91668
age_at_survey:kelp_annual_100                -0.02467    0.03663 1744.05339
hard_bottom_500:site_typeMPA:age_at_survey    0.34609    0.04690 1739.21871
site_typeMPA:age_at_survey:kelp_annual_100    0.01675    0.04669 1730.99960
                                           t value Pr(>|t|)    
(Intercept)                                  0.720 0.523638    
hard_bottom_500                             -1.683 0.092580 .  
site_typeMPA                                 6.868 9.01e-12 ***
age_at_survey                                0.883 0.377473    
kelp_annual_100                             -1.528 0.126611    
hard_bottom_500:site_typeMPA                 9.542  < 2e-16 ***
hard_bottom_500:age_at_survey                0.943 0.345589    
site_typeMPA:age_at_survey                   3.022 0.002545 ** 
site_typeMPA:kelp_annual_100                 3.722 0.000204 ***
age_at_survey:kelp_annual_100               -0.674 0.500710    
hard_bottom_500:site_typeMPA:age_at_survey   7.379 2.45e-13 ***
site_typeMPA:age_at_survey:kelp_annual_100   0.359 0.719863    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
             (Intr) hr__500 st_MPA ag_t_s k__100 hr__500:_MPA h__500:__
hrd_btt_500  -0.022                                                    
site_typMPA  -0.081 -0.016                                             
age_at_srvy   0.016  0.065  -0.055                                     
klp_nnl_100   0.021  0.122  -0.151  0.345                              
hr__500:_MPA -0.012 -0.304   0.012  0.040  0.014                       
hrd__500:__   0.007 -0.006  -0.022  0.123 -0.017  0.008                
st_tyMPA:__  -0.004 -0.077   0.025 -0.694 -0.237 -0.055       -0.085   
s_MPA:__100  -0.016 -0.125   0.114 -0.283 -0.768 -0.062        0.014   
ag_t_:__100   0.033  0.028  -0.220  0.276  0.577  0.027        0.065   
h__500:_MPA: -0.006  0.011  -0.037 -0.067  0.011  0.010       -0.481   
s_MPA:__:__  -0.024 -0.036   0.239 -0.213 -0.442 -0.064       -0.050   
             st_MPA:__ s_MPA:__1 a__:__ h__500:_MPA:
hrd_btt_500                                         
site_typMPA                                         
age_at_srvy                                         
klp_nnl_100                                         
hr__500:_MPA                                        
hrd__500:__                                         
st_tyMPA:__                                         
s_MPA:__100   0.269                                 
ag_t_:__100  -0.195    -0.451                       
h__500:_MPA: -0.038    -0.085    -0.033             
s_MPA:__:__   0.161     0.527    -0.761 -0.113      

Log+1 Transformed DV

log(B+1) ~ H*ST + K*ST + ST*A (No Scaling)

  • Removed year as it had zero variance (singular fit)
  • Significant: H*ST, K*ST, ST*A, ST
  • Almost significant: ST
  • Same conclusions as version without log-transformation (e.g. vs. glm1) except ST is almost significant now
glm5 <- lmer(log_kg_per_m2 ~ hard_bottom_500 * site_type + kelp_annual_100 * site_type + site_type * age_at_survey +
               (1|bioregion) + (1|affiliated_mpa), data = data_sp, REML = FALSE)
Warning: Some predictor variables are on very different scales: consider
rescaling
Warning: Some predictor variables are on very different scales: consider
rescaling
attr(glm5, "name") <- "GLMM: log(B+1) ~ H*ST + K*ST + ST*A (Not Scaled)"

summary(glm5)
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: log_kg_per_m2 ~ hard_bottom_500 * site_type + kelp_annual_100 *  
    site_type + site_type * age_at_survey + (1 | bioregion) +  
    (1 | affiliated_mpa)
   Data: data_sp

     AIC      BIC   logLik deviance df.resid 
-15575.4 -15515.2   7798.7 -15597.4     1737 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6738 -0.3367 -0.0312  0.1554 14.2151 

Random effects:
 Groups         Name        Variance  Std.Dev. 
 affiliated_mpa (Intercept) 9.194e-07 0.0009588
 bioregion      (Intercept) 1.949e-06 0.0013959
 Residual                   7.571e-06 0.0027516
Number of obs: 1748, groups:  affiliated_mpa, 23; bioregion, 3

Fixed effects:
                               Estimate Std. Error         df t value Pr(>|t|)
(Intercept)                   1.966e-03  9.007e-04  3.624e+00   2.182 0.101554
hard_bottom_500              -9.472e-10  5.746e-10  1.632e+03  -1.648 0.099466
site_typeMPA                 -3.112e-03  4.109e-04  1.721e+03  -7.575 5.82e-14
kelp_annual_100              -3.866e-08  2.891e-08  1.745e+03  -1.337 0.181362
age_at_survey                 1.915e-05  1.765e-05  1.743e+03   1.085 0.278031
hard_bottom_500:site_typeMPA  1.000e-08  1.068e-09  1.574e+03   9.362  < 2e-16
site_typeMPA:kelp_annual_100  1.516e-07  3.598e-08  1.745e+03   4.213 2.64e-05
site_typeMPA:age_at_survey    9.633e-05  2.538e-05  1.734e+03   3.795 0.000153
                                
(Intercept)                     
hard_bottom_500              .  
site_typeMPA                 ***
kelp_annual_100                 
age_at_survey                   
hard_bottom_500:site_typeMPA ***
site_typeMPA:kelp_annual_100 ***
site_typeMPA:age_at_survey   ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) hr__500 st_MPA k__100 ag_t_s h__500: s_MPA:__1
hrd_btt_500 -0.232                                               
site_typMPA -0.145  0.297                                        
klp_nnl_100 -0.124  0.130   0.197                                
age_at_srvy -0.184  0.061   0.361  0.246                         
h__500:_MPA  0.041 -0.308  -0.740 -0.003  0.033                  
s_MPA:__100  0.099 -0.124  -0.255 -0.732 -0.200 -0.037           
st_tyMPA:__  0.138 -0.074  -0.521 -0.161 -0.681 -0.050   0.215   
fit warnings:
Some predictor variables are on very different scales: consider rescaling

log(B+1) ~ H*ST + K*ST + ST*A (Scaled)

  • Year has zero variance (singular fit)
  • Significant: H*ST, K*ST, ST*A
  • Almost significant: ST
  • Same conclusions as version without log-transformation (e.g. vs. glm2) and same conclusions as non-scaled (e.g. vs. glm5)
glm6 <- lmer(log_kg_per_m2 ~ hard_bottom_500 * site_type + kelp_annual_100 * site_type + site_type * age_at_survey +
               (1|bioregion) + (1|affiliated_mpa), data = data_sp_scaled, REML = FALSE)

attr(glm6, "name") <- "GLMM: log(B+1) ~ H*ST + K*ST + ST*A (Scaled)"

summary(glm6)
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: log_kg_per_m2 ~ hard_bottom_500 * site_type + kelp_annual_100 *  
    site_type + site_type * age_at_survey + (1 | bioregion) +  
    (1 | affiliated_mpa)
   Data: data_sp_scaled

     AIC      BIC   logLik deviance df.resid 
  4301.3   4361.5  -2139.7   4279.3     1737 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6738 -0.3367 -0.0312  0.1554 14.2151 

Random effects:
 Groups         Name        Variance Std.Dev.
 affiliated_mpa (Intercept) 0.07978  0.2825  
 bioregion      (Intercept) 0.16910  0.4112  
 Residual                   0.65702  0.8106  
Number of obs: 1748, groups:  affiliated_mpa, 23; bioregion, 3

Fixed effects:
                               Estimate Std. Error         df t value Pr(>|t|)
(Intercept)                     0.19720    0.25377    3.04449   0.777 0.493008
hard_bottom_500                -0.04455    0.02703 1632.16947  -1.648 0.099466
site_typeMPA                    0.28889    0.04074 1742.21461   7.092 1.92e-12
kelp_annual_100                -0.04526    0.03385 1744.66293  -1.337 0.181362
age_at_survey                   0.03038    0.02800 1743.04739   1.085 0.278031
hard_bottom_500:site_typeMPA    0.47034    0.05024 1574.09919   9.362  < 2e-16
site_typeMPA:kelp_annual_100    0.17751    0.04213 1745.08538   4.213 2.64e-05
site_typeMPA:age_at_survey      0.15279    0.04026 1734.05702   3.795 0.000153
                                
(Intercept)                     
hard_bottom_500              .  
site_typeMPA                 ***
kelp_annual_100                 
age_at_survey                   
hard_bottom_500:site_typeMPA ***
site_typeMPA:kelp_annual_100 ***
site_typeMPA:age_at_survey   ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) hr__500 st_MPA k__100 ag_t_s h__500: s_MPA:__1
hrd_btt_500 -0.024                                               
site_typMPA -0.080 -0.007                                        
klp_nnl_100  0.003  0.130  -0.030                                
age_at_srvy  0.007  0.061   0.008  0.246                         
h__500:_MPA -0.013 -0.308   0.026 -0.003  0.033                  
s_MPA:__100 -0.003 -0.124  -0.021 -0.732 -0.200 -0.037           
st_tyMPA:__  0.003 -0.074  -0.024 -0.161 -0.681 -0.050   0.215   

log(B+1) ~ H*ST*A + K*ST*A (Not Scaled)

  • Removed year as it had zero variance (singular fit)
  • Significant: H*ST*A, K*ST, ST*A
  • Almost significant: H, ST
  • Same conclusions as version without log-transformation (e.g. vs. glm3)
glm7 <- lmer(log_kg_per_m2 ~ hard_bottom_500 * site_type * age_at_survey + kelp_annual_100 * site_type * age_at_survey + 
               (1|bioregion) + (1|affiliated_mpa), data = data_sp, REML = FALSE)
Warning: Some predictor variables are on very different scales: consider
rescaling
Warning: Some predictor variables are on very different scales: consider
rescaling
attr(glm7, "name") <- "GLMM: log(B+1) ~ H*ST*A + K*ST*A (Not Scaled)"

summary(glm7)
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: log_kg_per_m2 ~ hard_bottom_500 * site_type * age_at_survey +  
    kelp_annual_100 * site_type * age_at_survey + (1 | bioregion) +  
    (1 | affiliated_mpa)
   Data: data_sp

     AIC      BIC   logLik deviance df.resid 
-15651.1 -15569.1   7840.5 -15681.1     1733 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0105 -0.3230 -0.0405  0.1650 14.6102 

Random effects:
 Groups         Name        Variance  Std.Dev. 
 affiliated_mpa (Intercept) 9.051e-07 0.0009514
 bioregion      (Intercept) 2.093e-06 0.0014466
 Residual                   7.214e-06 0.0026859
Number of obs: 1748, groups:  affiliated_mpa, 23; bioregion, 3

Fixed effects:
                                             Estimate Std. Error         df
(Intercept)                                 2.144e-03  9.514e-04  3.941e+00
hard_bottom_500                            -1.661e-09  9.417e-10  1.744e+03
site_typeMPA                                4.850e-04  6.156e-04  1.739e+03
age_at_survey                              -2.221e-06  3.145e-05  1.740e+03
kelp_annual_100                            -2.039e-08  4.008e-08  1.746e+03
hard_bottom_500:site_typeMPA               -1.408e-09  1.854e-09  1.717e+03
hard_bottom_500:age_at_survey               8.557e-11  9.021e-11  1.736e+03
site_typeMPA:age_at_survey                 -3.384e-04  6.114e-05  1.737e+03
site_typeMPA:kelp_annual_100                1.332e-07  5.319e-08  1.737e+03
age_at_survey:kelp_annual_100              -3.939e-09  5.800e-09  1.744e+03
hard_bottom_500:site_typeMPA:age_at_survey  1.368e-09  1.849e-10  1.739e+03
site_typeMPA:age_at_survey:kelp_annual_100  2.633e-09  7.394e-09  1.731e+03
                                           t value Pr(>|t|)    
(Intercept)                                  2.254   0.0883 .  
hard_bottom_500                             -1.764   0.0778 .  
site_typeMPA                                 0.788   0.4308    
age_at_survey                               -0.071   0.9437    
kelp_annual_100                             -0.509   0.6110    
hard_bottom_500:site_typeMPA                -0.759   0.4479    
hard_bottom_500:age_at_survey                0.949   0.3430    
site_typeMPA:age_at_survey                  -5.535 3.58e-08 ***
site_typeMPA:kelp_annual_100                 2.504   0.0124 *  
age_at_survey:kelp_annual_100               -0.679   0.4972    
hard_bottom_500:site_typeMPA:age_at_survey   7.400 2.10e-13 ***
site_typeMPA:age_at_survey:kelp_annual_100   0.356   0.7218    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
             (Intr) hr__500 st_MPA ag_t_s k__100 hr__500:_MPA h__500:__
hrd_btt_500  -0.315                                                    
site_typMPA  -0.174  0.398                                             
age_at_srvy  -0.290  0.669   0.428                                     
klp_nnl_100  -0.130  0.118   0.185  0.309                              
hr__500:_MPA  0.118 -0.431  -0.889 -0.318 -0.053                       
hrd__500:__   0.233 -0.803  -0.353 -0.813 -0.094  0.405                
st_tyMPA:__   0.154 -0.350  -0.830 -0.501 -0.155  0.736        0.413   
s_MPA:__100   0.093 -0.088  -0.179 -0.229 -0.714 -0.040        0.068   
ag_t_:__100   0.054 -0.035  -0.108 -0.248 -0.708  0.043        0.065   
h__500:_MPA: -0.115  0.391   0.740  0.387  0.050 -0.826       -0.481   
s_MPA:__:__  -0.037  0.018   0.069  0.187  0.536  0.058       -0.050   
             st_MPA:__ s_MPA:__1 a__:__ h__500:_MPA:
hrd_btt_500                                         
site_typMPA                                         
age_at_srvy                                         
klp_nnl_100                                         
hr__500:_MPA                                        
hrd__500:__                                         
st_tyMPA:__                                         
s_MPA:__100   0.122                                 
ag_t_:__100   0.123     0.530                       
h__500:_MPA: -0.897     0.065    -0.033             
s_MPA:__:__  -0.060    -0.748    -0.761 -0.113      
fit warnings:
Some predictor variables are on very different scales: consider rescaling

log(B+1) ~ H*ST*A + K*ST*A (Scaled)

  • Removed year as it had zero variance (singular fit)
  • Significant: H*ST*A, H*ST, K*ST, ST*A, ST
  • Same conclusions as version without log-transformation (e.g. vs. glm4) but different from non-scaled version b/c ST now highly significant on its own, and H*ST is significant on its own, hard bottom on its own is nearly significant
glm8 <- lmer(log_kg_per_m2 ~ hard_bottom_500 * site_type * age_at_survey + kelp_annual_100 * site_type * age_at_survey + 
               (1|bioregion) + (1|affiliated_mpa), data = data_sp_scaled, REML = FALSE)

attr(glm8, "name") <- "GLMM: log(B+1) ~ H*ST*A + K*ST*A (Scaled)"

summary(glm8)
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: log_kg_per_m2 ~ hard_bottom_500 * site_type * age_at_survey +  
    kelp_annual_100 * site_type * age_at_survey + (1 | bioregion) +  
    (1 | affiliated_mpa)
   Data: data_sp_scaled

     AIC      BIC   logLik deviance df.resid 
  4225.6   4307.6  -2097.8   4195.6     1733 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0105 -0.3230 -0.0405  0.1650 14.6102 

Random effects:
 Groups         Name        Variance Std.Dev.
 affiliated_mpa (Intercept) 0.07854  0.2803  
 bioregion      (Intercept) 0.18158  0.4261  
 Residual                   0.62600  0.7912  
Number of obs: 1748, groups:  affiliated_mpa, 23; bioregion, 3

Fixed effects:
                                             Estimate Std. Error         df
(Intercept)                                   0.18905    0.26173    3.00929
hard_bottom_500                              -0.04457    0.02642 1637.48328
site_typeMPA                                  0.28265    0.04105 1741.42207
age_at_survey                                 0.02536    0.02861 1743.62208
kelp_annual_100                              -0.06233    0.04056 1743.89705
hard_bottom_500:site_typeMPA                  0.47054    0.04920 1587.09477
hard_bottom_500:age_at_survey                 0.02167    0.02285 1736.20266
site_typeMPA:age_at_survey                    0.12248    0.04035 1734.00634
site_typeMPA:kelp_annual_100                  0.18165    0.04868 1743.85438
age_at_survey:kelp_annual_100                -0.02483    0.03656 1744.00838
hard_bottom_500:site_typeMPA:age_at_survey    0.34648    0.04682 1739.16854
site_typeMPA:age_at_survey:kelp_annual_100    0.01660    0.04661 1730.96135
                                           t value Pr(>|t|)    
(Intercept)                                  0.722 0.522162    
hard_bottom_500                             -1.687 0.091773 .  
site_typeMPA                                 6.885 8.02e-12 ***
age_at_survey                                0.886 0.375603    
kelp_annual_100                             -1.537 0.124554    
hard_bottom_500:site_typeMPA                 9.564  < 2e-16 ***
hard_bottom_500:age_at_survey                0.949 0.342965    
site_typeMPA:age_at_survey                   3.035 0.002439 ** 
site_typeMPA:kelp_annual_100                 3.732 0.000196 ***
age_at_survey:kelp_annual_100               -0.679 0.497184    
hard_bottom_500:site_typeMPA:age_at_survey   7.400 2.10e-13 ***
site_typeMPA:age_at_survey:kelp_annual_100   0.356 0.721838    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
             (Intr) hr__500 st_MPA ag_t_s k__100 hr__500:_MPA h__500:__
hrd_btt_500  -0.022                                                    
site_typMPA  -0.080 -0.016                                             
age_at_srvy   0.016  0.065  -0.055                                     
klp_nnl_100   0.021  0.122  -0.151  0.345                              
hr__500:_MPA -0.012 -0.304   0.012  0.040  0.014                       
hrd__500:__   0.007 -0.006  -0.022  0.123 -0.017  0.008                
st_tyMPA:__  -0.004 -0.077   0.025 -0.694 -0.237 -0.055       -0.085   
s_MPA:__100  -0.016 -0.125   0.114 -0.283 -0.768 -0.063        0.014   
ag_t_:__100   0.033  0.028  -0.220  0.276  0.577  0.027        0.065   
h__500:_MPA: -0.006  0.010  -0.037 -0.067  0.011  0.010       -0.481   
s_MPA:__:__  -0.024 -0.036   0.239 -0.213 -0.442 -0.064       -0.050   
             st_MPA:__ s_MPA:__1 a__:__ h__500:_MPA:
hrd_btt_500                                         
site_typMPA                                         
age_at_srvy                                         
klp_nnl_100                                         
hr__500:_MPA                                        
hrd__500:__                                         
st_tyMPA:__                                         
s_MPA:__100   0.269                                 
ag_t_:__100  -0.195    -0.451                       
h__500:_MPA: -0.038    -0.085    -0.033             
s_MPA:__:__   0.161     0.527    -0.761 -0.113      

Log + Small Constant Transformation

log(B+c) ~ H*ST + K*ST + ST*A (No Scaling)

  • Now there is nonzero variance for year
  • H*ST and K*ST is still significant but ST*A is not anymore
  • Now intercept is significant, H alone is significant, Age alone is significant
  • Changes compared to the non-transformed and log+1 transformation!
glm9 <- lmer(log_kg_per_m2_c ~ hard_bottom_500 * site_type + kelp_annual_100 * site_type + site_type * age_at_survey +
               (1|bioregion) + (1|affiliated_mpa) + (1|year), data = data_sp, REML = FALSE)
Warning: Some predictor variables are on very different scales: consider
rescaling
Warning: Some predictor variables are on very different scales: consider
rescaling
attr(glm9, "name") <- "GLMM: log(B+c) ~ H*ST + K*ST + ST*A (Not Scaled)"

summary(glm9)
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: log_kg_per_m2_c ~ hard_bottom_500 * site_type + kelp_annual_100 *  
    site_type + site_type * age_at_survey + (1 | bioregion) +  
    (1 | affiliated_mpa) + (1 | year)
   Data: data_sp

     AIC      BIC   logLik deviance df.resid 
  9139.7   9205.3  -4557.9   9115.7     1736 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2706 -0.5521 -0.0645  0.5798  3.2244 

Random effects:
 Groups         Name        Variance Std.Dev.
 affiliated_mpa (Intercept)  2.3174  1.5223  
 year           (Intercept)  0.1354  0.3679  
 bioregion      (Intercept)  8.0802  2.8426  
 Residual                   10.2866  3.2073  
Number of obs: 1748, groups:  affiliated_mpa, 23; year, 21; bioregion, 3

Fixed effects:
                               Estimate Std. Error         df t value Pr(>|t|)
(Intercept)                  -1.213e+01  1.737e+00  3.365e+00  -6.986 0.004100
hard_bottom_500              -1.568e-06  6.769e-07  1.698e+03  -2.316 0.020661
site_typeMPA                 -1.091e+00  4.817e-01  1.723e+03  -2.265 0.023641
kelp_annual_100              -4.279e-05  3.403e-05  1.746e+03  -1.257 0.208803
age_at_survey                 1.008e-01  2.487e-02  5.776e+01   4.054 0.000152
hard_bottom_500:site_typeMPA  5.806e-06  1.259e-06  1.667e+03   4.610 4.33e-06
site_typeMPA:kelp_annual_100  1.386e-04  4.210e-05  1.730e+03   3.291 0.001017
site_typeMPA:age_at_survey   -2.798e-02  2.963e-02  1.716e+03  -0.944 0.345069
                                
(Intercept)                  ** 
hard_bottom_500              *  
site_typeMPA                 *  
kelp_annual_100                 
age_at_survey                ***
hard_bottom_500:site_typeMPA ***
site_typeMPA:kelp_annual_100 ** 
site_typeMPA:age_at_survey      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) hr__500 st_MPA k__100 ag_t_s h__500: s_MPA:__1
hrd_btt_500 -0.142                                               
site_typMPA -0.086  0.285                                        
klp_nnl_100 -0.075  0.127   0.192                                
age_at_srvy -0.111  0.048   0.298  0.203                         
h__500:_MPA  0.022 -0.291  -0.742  0.002  0.029                  
s_MPA:__100  0.060 -0.126  -0.251 -0.728 -0.166 -0.042           
st_tyMPA:__  0.084 -0.075  -0.516 -0.161 -0.566 -0.052   0.216   
fit warnings:
Some predictor variables are on very different scales: consider rescaling

log(B+c) ~ H*ST + K*ST + ST*A (Scaled)

  • Now there is nonzero variance for year
  • H*ST and K*ST is still significant but ST*A is not anymore
  • Now H alone is significant, Age alone is significant, ST is significant
  • Changes compared to the non-transformed and log+1 transformation, similar to the non-scaled log+c version (e.g. glmm #9) except intercept is not significant in this version
glm10 <- lmer(log_kg_per_m2_c ~ hard_bottom_500 * site_type + kelp_annual_100 * site_type + site_type * age_at_survey +
               (1|bioregion) + (1|affiliated_mpa) + (1|year), data = data_sp_scaled, REML = FALSE)

attr(glm10, "name") <- "GLMM: log(B+c) ~ H*ST + K*ST + ST*A (Scaled)"

summary(glm10)
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: log_kg_per_m2_c ~ hard_bottom_500 * site_type + kelp_annual_100 *  
    site_type + site_type * age_at_survey + (1 | bioregion) +  
    (1 | affiliated_mpa) + (1 | year)
   Data: data_sp_scaled

     AIC      BIC   logLik deviance df.resid 
  3935.5   4001.1  -1955.7   3911.5     1736 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2706 -0.5521 -0.0645  0.5798  3.2244 

Random effects:
 Groups         Name        Variance Std.Dev.
 affiliated_mpa (Intercept) 0.118046 0.34358 
 year           (Intercept) 0.006895 0.08304 
 bioregion      (Intercept) 0.411906 0.64180 
 Residual                   0.523924 0.72383 
Number of obs: 1748, groups:  affiliated_mpa, 23; year, 21; bioregion, 3

Fixed effects:
                               Estimate Std. Error         df t value Pr(>|t|)
(Intercept)                     0.39905    0.38581    3.15169   1.034 0.373731
hard_bottom_500                -0.05650    0.02439 1698.09694  -2.316 0.020659
site_typeMPA                    0.15380    0.03644 1718.90196   4.221 2.56e-05
kelp_annual_100                -0.03838    0.03053 1745.56360  -1.257 0.208789
age_at_survey                   0.12251    0.03022   57.75493   4.054 0.000153
hard_bottom_500:site_typeMPA    0.20922    0.04538 1666.64841   4.610 4.33e-06
site_typeMPA:kelp_annual_100    0.12428    0.03776 1729.73277   3.291 0.001017
site_typeMPA:age_at_survey     -0.03400    0.03600 1716.16888  -0.944 0.345069
                                
(Intercept)                     
hard_bottom_500              *  
site_typeMPA                 ***
kelp_annual_100                 
age_at_survey                ***
hard_bottom_500:site_typeMPA ***
site_typeMPA:kelp_annual_100 ** 
site_typeMPA:age_at_survey      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) hr__500 st_MPA k__100 ag_t_s h__500: s_MPA:__1
hrd_btt_500 -0.015                                               
site_typMPA -0.047 -0.007                                        
klp_nnl_100  0.001  0.127  -0.030                                
age_at_srvy  0.021  0.048   0.007  0.203                         
h__500:_MPA -0.008 -0.291   0.026  0.002  0.029                  
s_MPA:__100 -0.001 -0.126  -0.022 -0.728 -0.166 -0.042           
st_tyMPA:__  0.002 -0.075  -0.025 -0.161 -0.566 -0.052   0.216   

log(B+c) ~ H*ST*A + K*ST*A (Not Scaled)

  • Now there is nonzero variance for year
  • H*ST*A is still significant but less so; ST*A and K*ST is significant
  • Intercept, H alone are now both highly significant; ST no longer significant
  • Now significant: H*A
  • Changes compared to the non-transformed and log+1 transformation.
glm11 <- lmer(log_kg_per_m2_c ~ hard_bottom_500 * site_type * age_at_survey + kelp_annual_100 * site_type * age_at_survey + 
               (1|bioregion) + (1|affiliated_mpa) + (1|year), data = data_sp, REML = FALSE)
Warning: Some predictor variables are on very different scales: consider
rescaling
Warning: Some predictor variables are on very different scales: consider
rescaling
attr(glm11, "name") <- "GLMM: log(B+c) ~ H*ST*A + K*ST*A (Not Scaled)"

summary(glm11)
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: log_kg_per_m2_c ~ hard_bottom_500 * site_type * age_at_survey +  
    kelp_annual_100 * site_type * age_at_survey + (1 | bioregion) +  
    (1 | affiliated_mpa) + (1 | year)
   Data: data_sp

     AIC      BIC   logLik deviance df.resid 
  9129.2   9216.7  -4548.6   9097.2     1732 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1972 -0.5341 -0.0619  0.5650  3.3157 

Random effects:
 Groups         Name        Variance Std.Dev.
 affiliated_mpa (Intercept)  2.3012  1.5170  
 year           (Intercept)  0.1151  0.3392  
 bioregion      (Intercept)  8.1154  2.8487  
 Residual                   10.1863  3.1916  
Number of obs: 1748, groups:  affiliated_mpa, 23; year, 21; bioregion, 3

Fixed effects:
                                             Estimate Std. Error         df
(Intercept)                                -1.162e+01  1.759e+00  3.509e+00
hard_bottom_500                            -3.469e-06  1.129e-06  1.746e+03
site_typeMPA                                1.612e-01  7.352e-01  1.733e+03
age_at_survey                               4.065e-02  3.939e-02  3.986e+02
kelp_annual_100                            -3.321e-05  4.814e-05  1.733e+03
hard_bottom_500:site_typeMPA                1.732e-06  2.220e-06  1.732e+03
hard_bottom_500:age_at_survey               2.276e-07  1.080e-07  1.727e+03
site_typeMPA:age_at_survey                 -1.825e-01  7.280e-02  1.722e+03
site_typeMPA:kelp_annual_100                1.637e-04  6.335e-05  1.724e+03
age_at_survey:kelp_annual_100              -2.962e-06  6.952e-06  1.739e+03
hard_bottom_500:site_typeMPA:age_at_survey  4.974e-07  2.204e-07  1.730e+03
site_typeMPA:age_at_survey:kelp_annual_100 -3.620e-06  8.800e-06  1.718e+03
                                           t value Pr(>|t|)   
(Intercept)                                 -6.606  0.00427 **
hard_bottom_500                             -3.072  0.00216 **
site_typeMPA                                 0.219  0.82644   
age_at_survey                                1.032  0.30261   
kelp_annual_100                             -0.690  0.49031   
hard_bottom_500:site_typeMPA                 0.780  0.43548   
hard_bottom_500:age_at_survey                2.109  0.03513 * 
site_typeMPA:age_at_survey                  -2.507  0.01227 * 
site_typeMPA:kelp_annual_100                 2.583  0.00987 **
age_at_survey:kelp_annual_100               -0.426  0.67009   
hard_bottom_500:site_typeMPA:age_at_survey   2.257  0.02413 * 
site_typeMPA:age_at_survey:kelp_annual_100  -0.411  0.68083   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
             (Intr) hr__500 st_MPA ag_t_s k__100 hr__500:_MPA h__500:__
hrd_btt_500  -0.203                                                    
site_typMPA  -0.111  0.391                                             
age_at_srvy  -0.186  0.628   0.403                                     
klp_nnl_100  -0.084  0.115   0.181  0.296                              
hr__500:_MPA  0.074 -0.422  -0.890 -0.298 -0.050                       
hrd__500:__   0.149 -0.802  -0.351 -0.765 -0.093  0.402                
st_tyMPA:__   0.099 -0.347  -0.829 -0.473 -0.152  0.735        0.411   
s_MPA:__100   0.060 -0.087  -0.175 -0.219 -0.711 -0.044        0.067   
ag_t_:__100   0.035 -0.032  -0.107 -0.239 -0.709  0.042        0.063   
h__500:_MPA: -0.074  0.388   0.740  0.364  0.047 -0.825       -0.479   
s_MPA:__:__  -0.024  0.016   0.067  0.179  0.533  0.060       -0.048   
             st_MPA:__ s_MPA:__1 a__:__ h__500:_MPA:
hrd_btt_500                                         
site_typMPA                                         
age_at_srvy                                         
klp_nnl_100                                         
hr__500:_MPA                                        
hrd__500:__                                         
st_tyMPA:__                                         
s_MPA:__100   0.120                                 
ag_t_:__100   0.122     0.527                       
h__500:_MPA: -0.897     0.068    -0.032             
s_MPA:__:__  -0.058    -0.747    -0.757 -0.115      
fit warnings:
Some predictor variables are on very different scales: consider rescaling

log(B+c) ~ H*ST*A + K*ST*A (Scaled)

  • Now there is nonzero variance for year
  • Significant: H*ST*A (but less so), H*ST, K*ST; now: H, ST, A, H*A
  • Changes compared to the non-transformed and log+1 transformation. ST*A no longer significant, and H, ST, A, H*A are now significant on their own.
glm12 <- lmer(log_kg_per_m2_c ~ hard_bottom_500 * site_type * age_at_survey + kelp_annual_100 * site_type * age_at_survey + 
               (1|bioregion) + (1|affiliated_mpa) + (1|year), data = data_sp_scaled, REML = FALSE)

attr(glm12, "name") <- "GLMM: log(B+c) ~ H*ST*A + K*ST*A (Scaled)"

summary(glm12)
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: log_kg_per_m2_c ~ hard_bottom_500 * site_type * age_at_survey +  
    kelp_annual_100 * site_type * age_at_survey + (1 | bioregion) +  
    (1 | affiliated_mpa) + (1 | year)
   Data: data_sp_scaled

     AIC      BIC   logLik deviance df.resid 
  3925.0   4012.4  -1946.5   3893.0     1732 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1972 -0.5341 -0.0619  0.5650  3.3157 

Random effects:
 Groups         Name        Variance Std.Dev.
 affiliated_mpa (Intercept) 0.117206 0.34235 
 year           (Intercept) 0.005862 0.07656 
 bioregion      (Intercept) 0.413363 0.64293 
 Residual                   0.518815 0.72029 
Number of obs: 1748, groups:  affiliated_mpa, 23; year, 21; bioregion, 3

Fixed effects:
                                             Estimate Std. Error         df
(Intercept)                                   0.39457    0.38632    3.14695
hard_bottom_500                              -0.05657    0.02429 1698.00635
site_typeMPA                                  0.14844    0.03743 1717.60963
age_at_survey                                 0.12560    0.03047   66.06507
kelp_annual_100                              -0.05195    0.03720 1743.96625
hard_bottom_500:site_typeMPA                  0.21192    0.04527 1669.08203
hard_bottom_500:age_at_survey                 0.04417    0.02095 1726.98524
site_typeMPA:age_at_survey                   -0.04867    0.03677 1715.10015
site_typeMPA:kelp_annual_100                  0.11971    0.04443 1726.04778
age_at_survey:kelp_annual_100                -0.01431    0.03358 1739.06377
hard_bottom_500:site_typeMPA:age_at_survey    0.09651    0.04276 1729.97794
site_typeMPA:age_at_survey:kelp_annual_100   -0.01748    0.04250 1718.08058
                                           t value Pr(>|t|)    
(Intercept)                                  1.021 0.379075    
hard_bottom_500                             -2.329 0.019983 *  
site_typeMPA                                 3.966 7.60e-05 ***
age_at_survey                                4.122 0.000107 ***
kelp_annual_100                             -1.397 0.162692    
hard_bottom_500:site_typeMPA                 4.681 3.08e-06 ***
hard_bottom_500:age_at_survey                2.109 0.035127 *  
site_typeMPA:age_at_survey                  -1.323 0.185850    
site_typeMPA:kelp_annual_100                 2.694 0.007123 ** 
age_at_survey:kelp_annual_100               -0.426 0.670088    
hard_bottom_500:site_typeMPA:age_at_survey   2.257 0.024134 *  
site_typeMPA:age_at_survey:kelp_annual_100  -0.411 0.680836    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
             (Intr) hr__500 st_MPA ag_t_s k__100 hr__500:_MPA h__500:__
hrd_btt_500  -0.014                                                    
site_typMPA  -0.050 -0.016                                             
age_at_srvy   0.024  0.053  -0.046                                     
klp_nnl_100   0.013  0.121  -0.150  0.292                              
hr__500:_MPA -0.008 -0.288   0.012  0.036  0.018                       
hrd__500:__   0.005 -0.008  -0.021  0.124 -0.020  0.009                
st_tyMPA:__  -0.002 -0.078   0.024 -0.597 -0.236 -0.058       -0.086   
s_MPA:__100  -0.010 -0.126   0.112 -0.239 -0.764 -0.066        0.016   
ag_t_:__100   0.020  0.030  -0.219  0.233  0.575  0.029        0.063   
h__500:_MPA: -0.004  0.009  -0.037 -0.065  0.010  0.006       -0.479   
s_MPA:__:__  -0.014 -0.038   0.239 -0.179 -0.440 -0.063       -0.048   
             st_MPA:__ s_MPA:__1 a__:__ h__500:_MPA:
hrd_btt_500                                         
site_typMPA                                         
age_at_srvy                                         
klp_nnl_100                                         
hr__500:_MPA                                        
hrd__500:__                                         
st_tyMPA:__                                         
s_MPA:__100   0.269                                 
ag_t_:__100  -0.194    -0.448                       
h__500:_MPA: -0.036    -0.084    -0.032             
s_MPA:__:__   0.161     0.526    -0.757 -0.115      

Quick overview:

  • Not much difference between scaled and non-scaled versions in terms of conclusions
  • Not much difference between non-transformed and log+1 transformed versions in terms of conclusions
  • Clear differences in the log+c transformed versions compared to the other two.
  • Next: Need to see whether scaling or the transformations yields better fit.

Comparison & Diagnostics

Partial Residual Plots

B ~ H*ST + K*ST + ST*A (Not Scaled)

plot(predictorEffects(glm1, ~ age_at_survey, partial.residuals = TRUE),
     id = list(n = 1), residuals.pch = 19, residuals.cex = 0.2)

plot(predictorEffects(glm1, ~ hard_bottom_500, partial.residuals = TRUE), 
     id = list(n = 1), residuals.pch = 19, residuals.cex = 0.2)

plot(predictorEffects(glm1, ~ kelp_annual_100, partial.residuals = TRUE),
     axes = list(x = list(rotate = 25)), 
     id = list(n = 1), residuals.pch = 19, residuals.cex = 0.2)

# Error likely means that the extremes in the hard bottom value are causing challenges?
# Looks like large number of low values in kelp canopy cover .... interesting.

B ~ H*ST + K*ST + ST*A (Scaled)

plot(predictorEffects(glm2, ~ age_at_survey, partial.residuals = TRUE),
     id = list(n = 1), residuals.pch = 19, residuals.cex = 0.2)

plot(predictorEffects(glm2, ~ hard_bottom_500, partial.residuals = TRUE),
     id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

plot(predictorEffects(glm2, ~ kelp_annual_100, partial.residuals = TRUE), 
     id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

B ~ H*ST*A + K*ST*A (Not Scaled)

plot(predictorEffects(glm3, ~ age_at_survey, partial.residuals = TRUE),
     axes = list(x = list(cex = 0.8), y = list(cex = 0.8, lim = c(-0.001, 0.02))),
     lattice = list(strip = list(cex = 0.5)),
     id = list(n = 1, cex = 0.5), residuals.pch = 19, residuals.cex = 0.2)

plot(predictorEffects(glm3, ~ hard_bottom_500, partial.residuals = TRUE),
     axes = list(x = list(rotate = 25)), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

plot(predictorEffects(glm3, ~ kelp_annual_100, partial.residuals = TRUE),
     axes = list(x = list(rotate = 35)), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

B ~ H*ST*A + K*ST*A (Scaled)

plot(predictorEffects(glm4, ~ age_at_survey, partial.residuals = TRUE),
     axes = list( x = list(rotate = 25, cex = 0.8), y = list(cex = 0.8)),
     lattice = list(strip = list(cex = 0.5)),
     id = list(n = 1, cex = 0.5), residuals.pch = 19, residuals.cex = 0.2)

plot(predictorEffects(glm4, ~ hard_bottom_500, partial.residuals = TRUE),
     axes = list(x = list(rotate = 25), y = list(lim = c(-1, 5))), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

plot(predictorEffects(glm4, ~ kelp_annual_100, partial.residuals = TRUE),
     axes = list(x = list(rotate = 25), y = list(lim = c(-1, 5))), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

log(B+1) ~ H*ST + K*ST + ST*A (Not Scaled)

plot(predictorEffects(glm5, ~ age_at_survey, partial.residuals = TRUE),
     id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

plot(predictorEffects(glm5, ~ hard_bottom_500, partial.residuals = TRUE),
     axes = list(x = list(rotate = 25)), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

plot(predictorEffects(glm5, ~ kelp_annual_100, partial.residuals = TRUE),
     axes = list(x = list(rotate = 25)), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

# Error likely means that the extremes in the hard bottom value are causing challenges?

log(B+1) ~ H*ST + K*ST + ST*A (Scaled)

plot(predictorEffects(glm6, ~ age_at_survey, partial.residuals = TRUE), 
     id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

plot(predictorEffects(glm6, ~ hard_bottom_500, partial.residuals = TRUE), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

plot(predictorEffects(glm6, ~ kelp_annual_100, partial.residuals = TRUE), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

log(B+1) ~ H*ST*A + K*ST*A (Not Scaled)

plot(predictorEffects(glm7, ~ age_at_survey, partial.residuals = TRUE),
     axes = list( x = list(cex = 0.8), y = list(cex = 0.8)),
     lattice = list(strip = list(cex = 0.5)),
     id = list(n = 1, cex = 0.5), residuals.pch = 19, residuals.cex = 0.2)

plot(predictorEffects(glm7, ~ hard_bottom_500, partial.residuals = TRUE),
     axes = list(x = list(rotate = 25)), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

plot(predictorEffects(glm7, ~ kelp_annual_100, partial.residuals = TRUE),
     axes = list(x = list(rotate = 25)), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

log(B+1) ~ H*ST*A + K*ST*A (Scaled)

plot(predictorEffects(glm8, ~ age_at_survey, partial.residuals = TRUE),
     axes = list( x = list(rotate = 25, cex = 0.8), y = list(cex = 0.8)),
     lattice = list(strip = list(cex = 0.5)),
     id = list(n = 1, cex = 0.5), residuals.pch = 19, residuals.cex = 0.2)

plot(predictorEffects(glm8, ~ hard_bottom_500, partial.residuals = TRUE),
     axes = list(x = list(rotate = 25)), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

plot(predictorEffects(glm8, ~ kelp_annual_100, partial.residuals = TRUE),
     axes = list(x = list(rotate = 25)), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

log(B+c) ~ H*ST + K*ST + ST*A (Not Scaled)

plot(predictorEffects(glm9, ~ age_at_survey, partial.residuals = TRUE),
     axes = list(x = list(rotate = 25)), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

plot(predictorEffects(glm9, ~ hard_bottom_500, partial.residuals = TRUE),
     axes = list(x = list(rotate = 25)), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

plot(predictorEffects(glm9, ~ kelp_annual_100, partial.residuals = TRUE),
     axes = list(x = list(rotate = 25)), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

log(B+c) ~ H*ST + K*ST + ST*A (Scaled)

plot(predictorEffects(glm10, ~ age_at_survey, partial.residuals = TRUE),
     id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

plot(predictorEffects(glm10, ~ hard_bottom_500, partial.residuals = TRUE), 
     id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

plot(predictorEffects(glm10, ~ kelp_annual_100, partial.residuals = TRUE), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

log(B+c) ~ H*ST*A + K*ST*A (Not Scaled)

plot(predictorEffects(glm11, ~ age_at_survey, partial.residuals = TRUE),
     axes = list( x = list(cex = 0.8), y = list(cex = 0.8)),
     lattice = list(strip = list(cex = 0.5)),
     id = list(n = 1, cex = 0.5), residuals.pch = 19, residuals.cex = 0.2)

plot(predictorEffects(glm11, ~ hard_bottom_500, partial.residuals = TRUE),
     axes = list(x = list(rotate = 25)), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

plot(predictorEffects(glm11, ~ kelp_annual_100, partial.residuals = TRUE),
     axes = list(x = list(rotate = 25)), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

log(B+c) ~ H*ST*A + K*ST*A (Scaled)

plot(predictorEffects(glm12, ~ age_at_survey, partial.residuals = TRUE),
      axes = list( x = list(cex = 0.8), y = list(cex = 0.8)),
     lattice = list(strip = list(cex = 0.5)),
     id = list(n = 1, cex = 0.5), residuals.pch = 19, residuals.cex = 0.2)

plot(predictorEffects(glm12, ~ hard_bottom_500, partial.residuals = TRUE),
     axes = list(x = list(rotate = 25), y = list(lim = c(-1, 2))), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

plot(predictorEffects(glm12, ~ kelp_annual_100, partial.residuals = TRUE),
     axes = list(x = list(rotate = 25), y = list(lim = c(-1, 2))), id = list(n = 1),
     residuals.pch = 19,   
     residuals.cex = 0.2)

GLMM Summary

  • The diagnostic plots don’t show drastic differences among the non-transformed vs. log+1 transformed. Within those two, the non-scaled, 3-way interaction tends to have the lowest AIC.
  • The log+c transformation causes challenges because the zeroes become negative values (log of a value less than 1 is negative).
  • Scaling the kelp canopy cover may be creating challenges (scaled values go from -0.5 to 6… suggests lots of zeroes and some very high values). This could be because kelp is scaled across all years, rather than within each year?

Gamma GLMs

1. Fit Gamma GLM + 2way

# gglm1 <- glmer(kg_per_m2_c ~ hard_bottom_500 * site_type + kelp_annual_100 * site_type + site_type * age_at_survey + 
#                  (1|bioregion) + (1|affiliated_mpa) + (1|year),
#                family = Gamma(link = "log"), data = data_sp)

# does not fit